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Abstract 

The nature of freezing and melting transitions for a system of model colloids 
interacting by a DLVO potential in a spatially periodic external potential is 
studied using extensive Monte Carlo simulations. Detailed finite size scaling 
analyses of various thermodynamic quantities like the order parameter, its 
cumulants etc. are used to map the phase diagram of the system for various 
values of the reduced screening length KOg and the amplitude of the external 
potential. We find clear indication of a reentrant liquid phase over a significant 
region of the parameter space. Our simulations therefore show that the system 
of soft disks behaves in a fashion similar to charge stabilized colloids which are 
known to undergo an initial freezing, followed by a re-melting transition as the 
amplitude of the imposed, modulating field produced by crossed laser beams 
is steadily increased. Detailed analysis of our data shows several features 
consistent with a recent dislocation unbinding theory of laser induced melting. 
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I. INTRODUCTION 



The liquid-solid transition in two dimensional systems of particles under the influence of 
external modulating potentials has recently attracted a fair amount of attention from exper- 
iments [1-7], theory [8,9] and computer simulations [10-13]. This is partly due to the fact 
that well controlled, clean experiments can be performed using colloidal particles [16] con- 
fined between glass plates (producing essentially a two-dimensional system) and subjected 
to a spatially periodic electromagnetic field generated by interfering two, or more, crossed 
laser beams. One of the more surprising results of these studies, where a commensurate, 
one dimensional, modulating potential is imposed, is the fact that there exist regions in 
the phase diagram over which one observes reentrant [4-6] freezing/melting behavior. As 
a function of the laser field intensity the system first freezes from a modulated liquid to a 
two dimensional triangular solid. A further increase of the intensity confines the particles 
strongly within the troughs of the external potential, suppressing fiuctuations perpendicular 
to the troughs, which leads to an uncoupling of neighboring troughs and to re-melting. 

Our present understanding of this curious phenomenon has come from early mean field 
density functional [8] and more recent dislocation unbinding [9] calculations. The mean field 
theories neglect fiuctuations and therefore cannot explain reentrant behavior. The order 
of the transition is predicted to be first order for small laser field intensities, though for 
certain combinations of external potentials (which includes the specific geometry studied 
in the experiments and in this paper) the transition may become second order after going 
through a tricritical point. In general, though mean field theories are applicable in any 
dimension, the results are expected to be accurate only for higher dimensions and long 
ranged potentials. The vahdity of the predictions of such theories for the system under 
consideration is, therefore, in doubt. 

A more recent theory [9] extends the dislocation unbinding mechanism for two- 
dimensional melting [24,25] to systems under external potentials. For a two-dimensional 
triangular solid subjected to an external one-dimensional modulating potential, the only 
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dislocations involved are those which have their Burger's vectors parallel to the troughs of 
the potential. The system, therefore, maps onto an anisotropic, scalar Coulomb gas (or XY 
model) [9] in contrast to a vector Coulomb gas [24,25] for the pure 2D melting problem. 
Once bound dislocation pairs are integrated out, the melting temperature is obtained as 
a function of the renormalized or "effective" elastic constants which depend on external 
parameters like the strength of the potential, temperature and/or density. Though explicit 
calculations are possible only near the two extreme limits of zero and infinite field intensities, 
one can argue effectively that a reentrant melting transition is expected on general grounds 
quite independent of the detailed nature of the interaction potential for any two-dimensional 
system subject to such external potentials. The actual extent of this region could, of course, 
vary from system to system. In addition, these authors predict that the autocorrelation 
function of the Fourier components of the density (the Debye- Waller correlation function) 
decays algebraically in the solid phase at the transition with an universal exponent which 
depends only on the geometry and the magnitude of the reciprocal lattice vector. 

Computer simulation results in this field have so far been inconclusive. Early simula- 
tions [10] involving colloidal particles interacting via the Derjaguin, Landau, Verwey and 
Overbeek (DLVO) potential [16] found a large reentrant region in apparent agreement with 
later experiments. On closer scrutiny, though, quantitative agreement between simulation 
and experiments on the same system (but with slightly different parameters) appears to be 
poor [6]. Subsequent simulations [11-13] have questioned the findings of the earlier compu- 
tation and the calculated phase diagram does not show a significant reentrant liquid phase. 

Motivated, in part, by this controversy, in Ref. [14] we have recently investigated the 
freezing/melting behavior of a two dimensional hard disk system in an external potential. 
The pure hard disk system is rather well studied [17-20] by now and the nature of the 
melting transition in the absence of external potentials reasonably well explored. Also, 
there exist colloidal systems with hard interactions [16] so that, at least in principle, actual 
experiments using this system are possible. Finally, a hard disk simulation is relatively 
cheap to implement and one can make detailed studies of large systems without straining 
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computational resources. By these calculations we obtained a clear signature of a reentrant 
liquid phase showing that this phenomenon is indeed a general one as indicated in Ref. [9]. 

In the present paper we studied a system of particles interacting by a DLVO potential 
in a external periodic potential, motivated on one hand whether this reentrance scenario is 
dependent on the range of interaction, and on the other hand to compare it with experimental 
results [6]. 

The phase diagram has been computed by an application of finite size scaling methods 
similar to the methods used in our study of the hard disk systems in external potentials [14]. 

The rest of our paper is organized as follows. In Section II we specify the model and 
the simulation method including details of the finite size analysis used. In Section III we 
present our results for both zero and non-zero external potential, in particular results for the 
order parameter and its cumulants with a discussion on finite size effects. We also present 
other quantities like order parameter susceptibility, correlation functions and heat capacity, 
which further illustrates the nature of the phase transitions in this system. In Section IV 
we discuss our work in relation to the existing literature on this subject, summarize and 
conclude. 



II. MODEL AND METHOD 



A. The Model 



1. Potentials 

We study a system of N soft disks in a two dimensional box of fixed volume interacting 
with the DLVO pair potential 4>{fij) [16] between particles i and j with distance r^j. 



where R is the radius of the particles, = y eoe^ksT ^i^i^f inverse debye screening 

length, Z* is the effective surface charge, and is the dielectric constant of water. We used 



Sr — 78 and, unless otherwise indicated, 2R — 1.07//m and Z* — 7800. Additionally we 
chose a temperature of T = 293. ISX, and the particle density such that the particle spacing 
of an ideal lattice is Gg = 2.52578/xm. We then obtain different values for the reduced inverse 
screening length kQs by varying k as needed. In our simulation we set 2R to be the unit 
length. The potential in eq. (1) mainly depends on the value of kGs, so all features found 
for this system should be valid also for slightly different values of the other parameters 
mentioned above. 

In addition a particle with coordinates {x, y) is exposed to an external periodic potential 
of the form: 

= sin(27rx/do) (2) 

The constant do in Eq.(2) is chosen such that, for a density p = N/S^Sy, the modulation 
is commensurate to a triangular lattice of disks with nearest neighbor distance a^: do — 
a. a/3/2. 

The main parameters which define our system are kGs and the reduced potential strength 
Vo/kBT — Vq, where ks is the Boltzmann constant. 

2. box geometry 

All of the data (unless otherwise indicated) presented for Vq* < 0.2 are obtained by a 
simulation in a rectangular box of size Sx * Sy {Sx/Sy — and periodic boundary 

conditions in x- and y-direction, i.e. exactly as in [14]. We will refer to it as 'fixed box 
geometry' in the rest of the paper. 

For V^o* 7^ the external potential modulates the structure of the fluid and the particles 
form troughs oriented in the ^/-direction. In order to avoid unphysical results, for V^* > 0.2 
we mainly used a box with periodic boundary conditions in x-direction and movable walls 
in y-direction, see Fig. 1, and we will call this 'variable box geometry'. The simulation box 
volume is fixed as well as the side length Sx, but in x-direction the box is divided into slabs 
of width do, centered around the minima of the external periodic potential. The wall at the 
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end of each slab can move at most Ug upwards or downwards around its equilibrium position: 
|6| < tts, such that each slab has variable length between Sy — 2as and Sy + 2as in ^/-direction. 
The averaged box geometry still is Sx/ Sy = \/3/2 as for the fixed one. The constraint for 
neighboring walls is to have a distance less than as/2: |c| < as/2. The walls are hard, so no 
particle can cross them. This is indicated as thick sohd line in Fig. 1. To accommodate the 
particles in the box as well as possible, additional 'boundary' particles were placed in center 
(/ = do/2) behind each wall at a distance of e = as/2. The boundary particles interact 
with the particles in the box by the usual DLVO potential, but do not interact with each 
other. The motion of the walls is done with a Monte-Carlo procedure and keeps the volume 
constant, so we are still in the NVT (but variable shape) ensemble. 

The movable walls are chosen to give the system additional degrees of freedom to relax 
internal stress. We were lead to this geometry by some unphysical results when using the 
fixed box geometry and higher particle numbers {N > 4096). For a detailed discussion see 
end of section IIIB. 

B. The Method 

1. Numerical Details 

We perform NVT Monte Carlo (MC) simulations [21,22] for the system with interactions 
given by Eqs. (1) and (2) for various values of V^* and KOg. 

Averages < • > of observables have been obtained with the canonical measure. In order to 
obtain thermodynamic quantities for a range of system sizes, we analyzed various quantities 
within subsystems and used < • >l to denote averages in it. The subsystems are of size 
Lx X Ly where and Ly are chosen as Ly — Las and — Ly\/3/2 — Ldo consistent with 
the geometry of the triangular lattice. A subbox of size L = 3 as shown in Fig. 2 contains 
in average A^^ = = 9 particles. 

Most of the simulations described below have been done for a total system size of A" = 
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1024 and N — 4096 particles, additional ones with N — 400. Phase transitions have been 
studied in most cases by starting in the ordered solid and increasing kQs for fixed Vq. Runs 
where kGs is decreased were also performed for comparison. 

A typical simulation run with 10^ Monte Carlo steps (MCS) per particle (including 3 
xlO^ MCS for relaxation) took about 50 CPU hours on a PII/500 MHz PC. In addition 
to ordinary (local) MC moves we also used 'trough moves', by which particle placements 
in neighboring troughs are tried. Besides producing faster equilibration, including such 
moves ensures that at high V^* the formation of dislocations is not artificially hindered since 
particles can in effect bypass each other more easily — this is very unlikely with purely local 
MC moves. 



2. Observables 



The potential energy of the system per particle s is computed by: 



1 ^ 

hr.T ^ 



1=1 



^(j){rij) + V{xi,yi 



and the heat capacity per particle from the fluctuations of e* 



(3) 



^ AT < (£*- < e* >)2 > (4) 

Kb 

The nature of the fluid-solid phase transition in two dimensions has been a topic of con- 
troversy throughout the last forty years [18-20,23-27]. It is well known that true long range 
positional order is absent in the infinitely large system due to low energy long wavelength ex- 
citations so that translational correlations decay algebraically. According to the dislocation 
unbinding mechanism [24-26] the two dimensional solid (with quasi long ranged positional 
and long ranged orientational order) first melts into a "hexatic" phase with no positional 
order but with quasi long ranged orientational order signified by an algebraic decay of bond- 
orientational correlation. A second KT transition, driven by disclination unbinding, leads 
to melting of the hexatic into the liquid, where both the orientational and positional order 
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is short ranged. Therefore a useful order parameter in zero external field is the orientational 
order parameter. For a particle j located at fj we define the local orientational order: 



where Nh is the number of nearest neighbors, and 9ij the angle between the axis f; — fj and 
an arbitrary reference axis. For the total system we use: 



In an external periodic field given by Eq.(2), however, the bond orientational order 
parameter is nonzero even in the fluid phase [9,12]. This is because for Vq ^ we have 
now a "modulated" liquid, in which local hexagons consisting of the six nearest neighbors 
of a particle are automatically oriented by the external field. Thus < i/j^ > is nonzero both 
in the (modulated) liquid and the crystalline phase and it cannot be used to study phase 
transitions in this system. The order parameters corresponding to a solid phase are the 
Fourier components of the (non-uniform) density p(r) calculated at the reciprocal lattice 
points {G}. This (infinite) set of numbers are all zero (for G 7^ ) in an uniform liquid 
phase and nonzero in a solid. We restrict ourselves to the star consisting of the six smallest 
reciprocal lattice vectors of the two-dimensional triangular lattice. In the modulated liquid 
phase that is relevant to our system, the Fourier components corresponding to two out of 
these six vectors, viz. those in the direction perpendicular to the troughs of the external 
potential, are nonzero [8]. The other four components of this set consisting of the one in 
the direction Gi (as defined for the ideal crystal in Fig. 2), and those equivalent to it by 
symmetry, are zero in the (modulated) liquid and nonzero in the solid (if there is true long 
range order). We therefore use the following order parameter: 




1 ^ 



and as orientational correlation function: 
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where fj is the position vector of the j*'* particle. The corresponding susceptibihty xGi is: 

ksTxG, = [{{^PG,f) - (5) 

To measure the positional correlation, we chose the Debye- Waller correlation function which 
we define as follows: 

Cg^{R) = I < e'^i("(^)-"(°)) > I 

where R points to the elementary cell of the ideal lattice, and u{R) is the deviation of the 
actual particle position from the ideal lattice: r = R + u[R). In this case wc have chosen the 
direction of R to lie along the y axis (i.e. along the troughs of the potential). In the solid 
we expect this quantity do decay algebraically, i.e. C^{y) oc [9,24], where rj^ depends 

on the elastic constants. 

ipOi is sensitive to the phase transition where positional order is lost. Therefore, when 
decreasing V^* we expect the phase boundary to converge to the corresponding transition 
at zero external potential. But in contrast to Vq 7^ where the crystal is oriented by the 
external potential, at Iq* = it is only weakly fixed by the boundary conditions and can 
start to rotate, so we can not apply ipoi there. For this purpose we use a slightly modified 
positional order parameter tpci at Vq = 0: the phase information of i/jq (of course, before 
taking the absolute value) is used to determine the orientation of the crystal, and then 
a tilted coordinate system is used to compute ipoi- We applied the same method when 
calculating C^^ (y) at Vq = 0. 

We have determined phase transition points by the order parameter cumulant intersection 
method. The fourth order cumulant Ul oi the order parameter distribution is given by [28]: 

C'.W.-.) = 1-3^11^ (6) 

In order to distuingish between the cumulants of ipQ and ipoi , we denote them with U and 
Ul,G: respectively. In the hquid (short ranged order) C/^ ^ 1/3 and in the sohd (long range 
order) Ul ^ '2/3 for L ^ 00. In case of a continuous transition close to the transition point 
the cumulant is only a function of the ratio of the system size ~ Lug and the correlation 
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length ^: UL^Lcis/i)- Since ^ diverges at the critical point the cumulants for different system 
sizes intersect in one point: [/li(0) = t^LjlO) = ^* is a non-trivial value, i.e. C/* 7^ 1/3 
and U* 7^ 2/3. Even for first order transitions these cumulants intersect [29] though the 
value U* oi Ul at the intersection is not universal any more. The intersection point can, 
therefore, be taken as the phase boundary regardless of the order of the transition. This is 
useful since the order of the melting transition in 2D either in the absence [18-20,23-27] or 
with [8-13] external potentials is not unequivocally settled. And there is also another aspect 
in our special case: since the positional order correlation is predicted to decay algebraically 
in the solid phase (quasi long ranged order), the whole solid can be seen as consisting of a 
line or area of critical points with temperature- dependent critical indices r)^{T) [24]. In that 
case we expect the cumulants to merge at a nontrivial value at the onset of the solid phase 
instead of intersecting, yielding a line of intersections. We indeed observed this behavior 
for hard disks in high external potentials [14]. For the same reasons the same behavior is 
expected for the ■06-cumulants at the liquid-hexatic transition and in the hexatic phase [18]. 
Also the very similar 2d-XY-spin model shows this behavior when using the magnetization 
as order parameter [30]. 

Note that though the order parameter < V'Gi > decays to zero with increasing system 
size even in the 2-d solid (assuming quasi long ranged order there), the cumulants will stay 
at the non-trivial value regardless of L. So for L — > 00 there should be a jump from 1/3 to 
this nontrivial value when crossing the phase boundary from liquid to solid, which underlines 
the usefulness of cumulants. 

In order to map the phase diagram we systematically vary the system parameters V^* 
and KQs to detect order parameter cumulant intersection- or merging points which are then 
identified with the phase boundary. 
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III. RESULTS AND DISCUSSION 



A. Zero external potential 

In this section we analyse the system properties for zero external field. In particular we 
present results for the order parameter, the cumulants, the correlation functions and the 
heat capacity for different values of kUs- In these studies we used the fixed box geometry. 

In Fig. 3 the cumulant of the ■06 order parameter versus kQs is shown for different 
subsystem sizes. We identify the phase transition value of kqQs at about 14.42 by locating the 
cumulant intersection point. Since the positional order is not well defined in two-dimensional 
systems, the positional order parameter ipQ^ shows a strong system-size dependency, see 
Fig. 4. The cumulants of ipci intersect at a value of Kcds ~ 14.25, which is slightly smaller 
than KqQs- This is in agreement with a KTHNY two stage melting scenario, in which the solid 
and the fiuid phase are separated by a "hexatic" region in the phase diagram, in which the 
positional order is short ranged and the bond-orientational order is long ranged. Both effects 
are detected by the two order parameters, ipQ being sensitive for the bond-orientational order 
and ipGi on the positional order. Surprisingly, however, though in the case of the hexatic 
phase one expects the ■^e-cumulants to coincide, they obviously don't in Fig. 3 (see also 
[18]). 

The Debye- Waller correlation functions Cciiv) for different values of nag are shown in 
Fig. 5. At the transition value Kcds = 14.25 for = 4096 we find a power law dependency 
of Coiiy) from y with an exponent rjoi ~ 0.28 which is well within the predicted range of 
[1/4,1/3]. In Fig. 6 the orientational correlation function versus distance is shown. This 
function reveals a power law dependency of the bond-orientation correlations at KqQs, the 
exponent value is about 1/4. The value of the exponent rje at the transition has been 
predicted by the KTHNY theory to be 1/4 [9], which is in agreement with our results. 

In Fig. 7 (left) wc present the heat capacity data versus KQg for different system sizes. 
It is obvious that the heat capacity does not show a singularity as would be expected in 
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case of a first order or a conventional second-order transition. The peak maxima are not 
very sharp, but are roughly located close to the value of KqQs, where the ipQ-order parameter 
cumulants intersect. The peak maxima thus do not agree with the cumulant intersection 
point of the ■0G^-order parameter, which is again in agreement with the KTHNY scenario. 
We note that the identification of the phase transition point by the heat capacity maxima 
may result in misleading results on the location of the transition points in the phase diagram. 
In particular, for a smaller system of N = 400 particles we find that configurations wherein 
the crystal is rotated by a tilt angle of a = ±30° (a extracted from the phase information 
of ipe) may be present. Since this is incompatible with the box geometry it leads to a higher 
energy and a lower measured ipoi ■ The value of ipoi , on the other hand, is not appreciably 
altered. This is shown in the time evolution of that system in Fig. 8. We also show the 
configuration of a tilted crystal with a = 30° (3.000.000 MCS) in Fig. 9 (left), and of an 
'correctly' aligned crystal {a = 0°, 7.700.000 MCS) (right). 

B. External periodic potential 

In this section we analyse the system properties in the presence of a periodic external 
potential. The studies in this section are mainly done with variable box geometry. Compar- 
ative studies with fixed box geometry show that the new method does not lead to artificial 
features, but rather gives improved results. For more details see the discussion of the Debye- 
Waller correlation function at the end of this section. 

Two examples for the Acas-dependency of the V'Gi -order parameter and the cumulants for 
an external potential amplitude of Vq = 2 and Vq = 1000 are shown in Figs. 10 and 11. We 
note that the cumulant intersection, which can be clearly identified for Vq = 2 in Fig. 10 
is developing towards an intersection "hne" for Vq — 1000, a behavior which was found in 
case of the hard disk system in external potentials [14] as well as in related systems with a 
KT transition like the XY-spin-model [30]. Another example is shown in Fig. 12. There kQs 
is kept fixed at 15.3 and Vq is varied. The starting point at Vq — 0.2 is in the mod. hquid 
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phase, crosses slightly the solid ('laser induced freezing'), and re-enters the mod. hquid at 
higher Vq* ('reentrance'). This is already a first sign of a 'reentrant' phase transition scenario. 
For the phase diagram these cumulant intersection- or merging values were used. 

In Fig. 13 the cumulant intersection values are shown as a function of Vq for fixed box 
geometry. We observe that U* is not an universal number but, nevertheless, goes to a 
limiting value for large Vq [30]. 

The amount of hysteresis effects on the location of the transition point has been analyzed 
for the case of V^* = 2 by a time consuming reverse density quench simulation, in which a 
path in phase space was chosen in the direction opposite to the standard path. The results 
of this study are shown in Fig. 14. Comparing these results with the ones of Fig. 10 reveals 
quite close agreements showing that only small hysterese effects are present in the system. 

The ipoi order parameter susceptibilities Xd are shown in Fig. 15 versus for different 
system sizes. We note that close to the transition a maximum develops, the value of the 
maximum increasing with the system size. In Fig. 16 the susceptibilities of the largest 
subsystems (L = 32) as functions of Kag are compared for different values of Vq. Clearly the 
peak position for Vq = 2 is shifted to larger values compared to the cases with Vq — 0.2 and 
Vq — 1000. This feature is another sign of a "reentrant" phase transition scenario in the 
phase diagram. Compared to the cumulant intersection values, xgi maxima are located at 
slightly higher KUg. This may be due to finite size effects, which often show the feature that 
phase transition points in finite systems are shifted to slightly different values depending on 
the observable under investigation. In particular one expects (and we get) a shift towards 
parameter values in the disordered region (here a liquid, i.e. higher kQs) for the average 
order parameter and the susceptibility. 

In Fig. 7 (right) the heat capacity for Vq = 2 and different system sizes is shown. The 
peak is nearly independent from system size and shifted towards the liquid phase with respect 
to the order parameter cumulant intersection value, i.e. the same behavior as for V^* = 0. 

The advantage of the variable box geometry, especially for large systems, can be seen 
best by looking at the Debye- Waller correlation function. In Fig. 17 (left) an example is 
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shown for fixed box geometry. The crossing from the sohd phase with an algebraic decay to 
the mod. hquid with exponential decay is not monotonic, but at KUg — 15.7, Cciijj) drops 
to zero at 7/ = which is not physically meaningful. At a higher value = 16 it rises, 

and then falls again at nag — 16.4, showing a exponential decay as expected. In variable box 
geometry this feature doesn't show up, see Fig. 17 (right). Here we have a smooth transition 
from solid- to liquid- like behavior. We explain this strange behavior at nas = 15.7 in fixed 
geometry as follows: consider a system with N — 10000 particles. Without dislocations 
there will be an ideal lattice with Nf — 100 particles in each of the 100 troughs. Assuming 
dislocation unbinding as melting mechanism, consider the existence of some dislocations in 

— * 

the system with opposite burgers vectors b — ±056^^. One of these dislocations increases 
the number of particles in the troughs by one, while the other decreases it by one. We can 
have for example a situation where half of the troughs has 100 particles, and the other half 
has either 101 or 99 particles. If we now for simplicity assume that the distance between 
particles in a row is a = Ly/Nf, calculating the pair correlation function g{y) along a trough 
will yield two peaks around y = Ly/2: one centered exactly at y = Ly/2 from the troughs 
with Nt — 100 (a = a^), and the other centered at y = Ly/2 + as/2 due to the troughs with 
Nt — 99 or Nt — 101 (a 7^ a^). As consequence, Cciiv — Ly/2) will be zero. The same 
situation in the movable-walls geometry will not show these problems: the troughs with 101 
particles can expand a bit, while the ones with 99 particles can contract. Now in every trough 
is a = Qg, and Cci{y — Ly/2) is not necessarily zero. Also, the formation of a dislocation 
pair costs less energy and is closer to the true infinite system value. The discussion above 
is in some sense similar to the 2d-XY-spin model with a vortex in the center: in an infinite 
sample, the spins to the left and to the right have opposite spin directions, but periodic 
boundary conditions in a finite system will try to align them, so that the formation of the 
vortex is disturbed. Free boundary conditions won't cause this problem. 

In zero external potential with fixed geometry the formation of a dislocation pair is not 
so problematic, since the particles are not forced into troughs and have more degrees of 
freedom in movement. In the above example one end of a line of Nt — 101 particles could 
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make a slight shift in x-direction to access more space in y-direction. 

However, at the first data point in the sohd closest to the transition we find an algebraic 
decay Cgi (y) oc 1 /y*"?! ^ith rio, in the range of 0.25 . . . 0.34 for V^; = . . . 1000 and N = 1024 
particles. In [9] rja^ is predicted to be 1/4 at the transition. 

C. The Phase Diagram 

For each and Vq value we computed cumulants Ul^q for a range of subsystem sizes L 
and located intersection- or merging points which we identify with the phase boundary. We 
have obtained a detailed phase diagram for = 1024 particles which is shown in Fig. 18 
for fixed box geometry (left) and for variable box geometry (right). We want to emphasize 
that there are only slight differences and the general shape of the phase diagram is the same 
for both box geometries. At Vg* = ^-^so the V'e cumulant intersection value is plotted for 
comparison. The values of nag at the transition initially rise and subsequently drop as Vq 
increases. The maximum Kaas values are found for V^* f=:i 1 — 2. These transition points 
separate a high density solid from a low density modulated liquid. Thus, at a properly 
chosen kGs, we observe an initial freezing transition followed by a reentrant melting at a 
higher V^* value. Such an effect had been found earlier in experiments on colloidal systems 
in an external laser field [4-6] . 

In order to quantify residual finite size effects on the phase diagram, we have computed 
the transition points for different total system sizes. The resulting phase diagrams are 
shown in Fig. 19, again for fixed (left) and variable box geometry (right). We note that 
with increasing system size in fixed box geometry all transition points are slightly shifted 
to the solid region, whereas for the variable box this shift is towards the liquid region. One 
can see that the shift for the fixed box is much smaller at low Vq, and for the variable box 
it is much smaller at medium and high Vq". We also found that the cumulant intersection 
point smears out strongly if using the fixed box, N > 4096 and higher Vq, probably for the 
same reasons as those mentioned in the discussion of the Debye- Waller correlation function. 
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In the variable box there was no such problem. These features were the reason for us to 
use mainly the fixed box for < 0.2 and the variable box for Vq > 0.2. By the way, the 
same 'Debye- Waller problem' also occurs when simulating hard disks in external periodic 
potentials, N > 4096 and Vq medium or high, and can be solved again by using the variable 
box geometry. 

However, for all system sizes the structure of the phase diagram with a pronounced 
minimum at intermediate values of V^* is not affected by the shifts. 

The difference in the value of kqcis at the transition between the infinite and zero external 
potential cases, we find noasiVo — ~ '^GdsiVo = 0) ~ 0.82. This is not far away from 
0.608 which is a value predicted by [9]. 

We have also done simulations with slightly altered parameters, i.e using particles with 
diameter 2R = Siim, effective surface charge Z* = 20000 and = 8/im, to match the 
experiments in [6]. As expected, we only observe a slightly shift of the phase diagram of 
A{Kas) ~ 0.35 towards higher values of KUg. The experimental phase diagram [6] qualita- 
tively has the same shape as our results, but shows larger freezing- and reentrance regions 
and is shifted to higher values of nag at about A(«;as) f=:i 4.5 in average. The reasons for 
these differences are probably due to the particle interaction. We only use pairwise interac- 
tion, which is a good approximation for low particle densities [15]. But for higher particle 
densities many body interactions play a role because of macroion screening, which results 
in an effective pair potential that has considerable deviations [15] from a pure Yukawa-like 
potential like ours. In particular, there could be an attractive part. 

D. Scaling behavior 

We next try to determine the order of the phase transitions encountered in this system 
for two values of Vq. In order to investigate this issue we studied the scaling behavior of the 
order parameter, susceptibility and the order parameter cumulant near the phase boundary 
for a smaU (2) and a large (1000) Vq*. 
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Prom finite size scaling theory (for an overview see Ref. [22]) we expect these quantities 
to scale as [31] : 

< V^Gi >L L' ^ f{L/i) (7) 

XLhuTL-'^ ~ g{L/i) (8) 

Ul ~ h{L/^) (9) 

Here h = j3/u, c = ■y/u (for critical scaling), and /, g, h are scaling functions. Defining 
K — {Ktts — KGds)/ i^GCis-i wc cxpcct the Correlation length ^ to diverge as ^ oc for 
an ordinary critical point, while for a KT-transition we have an essential singularity and 
^ oc exp(aK~'^) when approaching the transition from the liquid side. 

In Fig. 20 we have plotted the left hand sides of Eqs.(7), (8) and (9) versus L/^ for 
Vq = 1000, where data points of the variable box geometry for 15.2 < nag < 16.0 have been 
considered and Koas = 15.1, obtained by cumulant intersection. In order not to introduce 
an unwarranted bias, we have separately considered ordinary critical scaling (left column) 
and a KT scaling form (right column) and adjusted the values of the parameters b, c and 
u, or a, b, c and i>, till we obtained collapse of our data onto a single curve determined by a 
least square estimator. ^ Good collapse of our data is observed for both scahng forms, the 
numerical values for D, 2b = r) and c = 2 — 77 for KT scaling (26 ^ 0.28, c fa 1.70, i> 0.37) 
are relatively close to the predicted values [9] (26 = rj = l/i, c = 1.75, i> = 0.5). The 
situation is similar for small Vq — 2, the quality of the collapse comparable to the one of 
Vq — 1000. The critical parameters were obtained in this case for values 15.4 < Kttg < 16.2, 
with Katts = 15.37. 

We have a good internal consistency between 77 = 0.25 . . . 0.33 extracted from the Debye- 
Waller correlation function, and the values obtained from data collapsing: r] — 0.28 for 



One remark concerning the KT-scaling: the errors in 1/ and a are relatively big because for 
example increasing i) and decreasing a by an appropriate amount resulted in a nearly equally good 
collapse. 
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Vq — 1000, and rj — 0.36 for Vq — 2. Our results for the numerical values of the parameters 
are summarized in Table I. 

A more precise classification of the phase transitions with the present data and system 
sizes is not easy. This topic is left for future work, in particular we plan to compute the elastic 
properties of the system by a method recently developed for the hard disk system [20,27] 
and to test the KT predictions [9]. 

IV. SUMMARY AND CONCLUSION 

In summary, we have calculated the phase diagram of a two dimensional system of soft 
disks, interacting via a DLVO potential, in an external sinusoidal potential. We find freezing 
followed by reentrant melting transitions over a significant region of the phase diagram in 
tune with results on hard disks [14], previous experiments on colloids [4-6] and with the 
expectations of a dislocation unbinding theory [9] . One of the main features of our calculation 
is the method used to locate phase boundaries. In contrast to earher simulations [10-13] 
which used either the jump of the order parameter or specific heat maxima to locate the 
phase transition, we used the more reliable cumulant intersection method. It must be noted 
that the specific heat in this system does not show a strong peak at the phase transition 
density so that its use may lead to confusing results. This, in our opinion, may be the 
reason for part of the controversy in this field. It is possible that earlier simulations which 
used smaller systems and no systematic finite size analysis may have overlooked this feature 
of cv which becomes apparent only in computations involving large system sizes. We have 
shown that finite size scaling of the order parameter cumulants as obtained from subsystem 
or subblock analysis, on the other hand, yields an accurate phase diagram. 

What is the order of the phase transition? We know that [18-20] for the pure hard 
disk system in two dimensions this question is quite difficult to answer and our present 
understanding [20] is that this system shows a KTHNY transition. In our system of soft 
disks, for zero external potential we can rule out a strong first order transition, although 
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smaller systems show a feature (double peak in the internal energy) which mimics such a 
behavior. We find several features which are consistent with the KT theory, but also one 
which is not. Upon turning on the external periodic potential, the difference between hexatic 
and liquid disappears, and an (anisotropic) KT transition [9] from the modulated liquid into 
the solid is expected. Our results show several features which suggest that this is what we 
have, but there are still some (not large) deviations from theory. Though we have discussed 
these observations in the rest of the paper, we list the important ones below for clarity: 

• The behavior of the cumulants near the transition is similar to an earlier work [30] on 
the XY system which shows a KT transition. 

• The specific heat is relatively featureless and does not scale with system size in a 
fashion expected of a true first order or conventional continuous transition. 

• The decay of the correlation functions is similar to what is predicted [9] for an 
anisotropic scalar Coulomb gas. 

• For two test values of Vq, the scaling of the order parameter, the susceptibility and 
the cumulant may be reasonable described by the KT theory. 

Of course, in order to resolve this issue unambiguously yet larger simulations are required. 
Also, we need to compute elastic properties [27,20] of this system in order to compare directly 
with the results of Ref. [9]. Work along these lines is in progress. 
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Table I Parameters in the scahng plots (see Fig. (20)) for Vq* = 2 and = 1000. The 

first three parameter columns are for critical scaling, the last four for KT scaling. The last 
line shows the predictions of KT theory. 
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FIGURES 
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FIG. 1. Schematic picture of the simulation box geometry used mainly for V^* > 0.2. 




FIG. 2. Schematic picture of the system geometry showing the direction Gi along which 
crystalline order develops at the transition modulated liquid to solid. The four vectors obtained 
by rotating Gi anti-clockwise by 60° and/or reflecting about the origin are equivalent. 
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FIG. 3. Cumulant of the tpQ order parameter versus KUg for various values of the system size L 
(N=4096, V^=0). 




FIG. 4. Average (left) and cumulant (right) of the ipOi order parameter versus KUg for various 

values of the system size L (N=4096, V^*=0). 
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FIG. 5. Debye Waller correlation function versus y for various values of KUg (N=4096, Vq = 0). 
Dashed line: Schematic picture of the functional decay with exponent 1/4, dotted line: schematic 
picture of the functional decay with exponent 1/3. 
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FIG. 6. Orientational correlation function versus distance for various values of kQs (N=4096, 
V^* = 0). Dashed line: schematic picture of the functional decay with exponent 1/4. 
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FIG. 7. Heat capacity versus for various numbers of particles. Left figure: V^* = computed 
with fixed box geometry, right figure: Vq = 2, computed with variable box geometry. 
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FIG. 8. System evolution versus Monte Carlo steps (N=400, V'o*=0, Ka5=14.4). Prom bottom 
to top: energy e*, jpQ order parameter, angle a of lattice direction, jpc order parameter, i/jq- 
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FIG. 9. Configurations after 3.000.000 MCS (left) and 7.700.000 MCS (right) from the system 
of Fig. 8. 
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FIG. 13. Cumulant intersection values of the V'Gi -order parameter versus V^* for fixed box 
geometry (N=1024). The data for variable box geometry is similar. 
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FIG. 14. Cumulant of the V'Gi order parameter versus nag for Vq = 2 and various system sizes 
L for a density quench path (N=4096). 
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FIG. 15. Susceptibility of the -order parameter versus KUg for Vq=2 and various values of 
the system size L (N=1024). 
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FIG. 16. Susceptibility of the -^g^ -order parameter versus kQs for various values of Vq (L=32, 
N=1024). 




y y 
FIG. 17. Debye- Waller correlation function Coiiy) versus y for V^* = 2 and various values 
of KQs (N=4096). Correlations for computations with fixed box geometry (left) and variable box 
geometry (right). 
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FIG. 18. Phase diagram. The points show the parameters for the cumulant intersection of the 
V'Gi -order parameter (N=1024). Left picture: Computations with fixed box geometry for all V^*. 
Right picture: computations with variable box geometry for all V^*. At V^* = the parameters for 
the cumulant intersection of the ipQ order parameter are shown for comparison. 
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FIG. 20. Scaling plots for the order parameter (first line), the order parameter cumulant (second 
line) and the order parameter susceptibility (third line) for Vq = 1000 assuming critical (left 
column) and KT scaling (right column). The total system size is N = 1024, for ^ we have used 
the expressions given after Eq.(9). 
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